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1.  Summary 

The  objective  of  this  research  is  to  develop  joint-inversion  methods  involving  P  travel 
times,  receiver  functions,  and  surface  wave  dispersion  measurements  and  to  apply  them 
to  the  western  China  region  to  obtain  3D  models  of  P  and  S  structures  of  the  crust  and 
upper  mantle.  Our  joint  inversion  methods  involve  obtaining  Pn  station  delays,  which  are 
then  used  with  surface  wave  dispersion  and  receiver  functions  at  each  station  to  invert 
jointly  for  S  structure,  the  Moho  depth,  and  Vp/Vs  ratios  in  the  crust.  Our  main 
accomplishments  include  the  following  four  areas.  (1)  We  successfully  implemented  a 
search-based  algorithm  (neighborhood  algorithm)  for  joint  inversion  of  surface  wave 
dispersion  data,  receiver  function,  and  Pn  delay  time.  The  implementation  uses  parallel 
programming  with  MPI  calls,  making  massive  data  processing  possible.  (2)  We  have 
systematically  carried  out  individual  components  of  the  project,  obtaining  unprecedented 
data  sets  for  surface  wave  dispersion,  receiver  functions,  and  Pn  tomography,  needed  for 
the  joint  inversion.  (3)  We  have  tested  the  joint  inversion  method  and  proposed  a 
practical  strategy  for  the  joint  inversions  of  the  real  data.  We  obtained  detailed  joint 
inversion  results  for  the  dense  HiClimb  array,  which  are  generally  consistent  with 
previous  results  but  show  considerable  difference  in  details  (crustal  structure  and  Moho 
depths).  (4)  We  have  performed  systematically  joint  inversions  of  all  available  stations  in 
western  China  and  obtained  a  3D  S  velocity  model,  an  average  crustal  Vp/Vs  map,  and 
derived  a  3D  P  velocity  model  of  the  region. 


2.  Introduction 

The  objective  of  this  research  is  to  develop  joint-inversion  methods  involving  P  travel 
times,  receiver  functions,  and  surface-wave  dispersion  measurements  from  ambient  noise 
correlation  and  to  apply  them  to  the  western  China  region  to  obtain  3D  models  of  P  and  S 
structures  of  the  crust  and  upper  mantle. 

In  seismic  inversion,  trade-offs  between  model  parameters  are  well  known.  To  resolve 
the  ambiguity  and  to  improve  resolution,  a  combination  of  different  data  sets  that  have 
sensitivities  to  different  parameters  is  required  or  a  priori  constraints  have  to  be  imposed. 
For  example,  it  has  long  been  recognized  that  teleseismic  P  receiver  functions  are 
primarily  sensitive  to  shear- wave  velocity  (Vs)  contrasts  and  the  depth- velocity  product, 
and  not  the  velocity  alone  (Ammon  et  al.,  1990).  On  the  other  hand,  surface  waves  are 
primarily  sensitive  to  the  vertical  shear-velocity  averages.  Thus,  joint  inversion  (or 
analysis)  of  receiver  functions  and  surface  wave  dispersion  has  now  been  commonly  used 
to  resolve  depth  resolution  of  Vs  structure  (e.g.,  Julia  et  al.  2000;  Ozalaybey  et  al.,  1997). 
Similarly,  crustal  thickness  estimated  using  only  the  delay  time  of  the  Moho  P-to-S 
converted  phase  (Ps)  relative  to  P  at  the  receiver  trades  off  strongly  with  the  crustal 
Vp/V s  ratio,  which  can  be  remedied  using  later  converted  phases  (Zhu  and  Kanamori, 
2000). 
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In  this  project,  we  obtain  and  exploit  a  vast  collection  of  data  that  we  have  access  to  for 
western  China  and  have  accumulated  over  decades.  The  data  includes  P  travel  times,  P 
receiver  functions,  Rayleigh  wave  group  and  phase  velocities  from  both  earthquakes  and 
ambient  noise.  We  then  combine  individual  data  components  into  a  joint  inversion  to 
derive  P  and  S  crustal  models.  In  particular,  we  add  additional  P  information  (Pn  delay 
times)  in  the  joint  inversion  of  receiver  functions  and  surface  wave  dispersion  to  obtain 
an  S  model  profile  and  average  crustal  Vp/Vs  ratios  for  stations  that  have  the  three  types 
of  data.  Essentially,  the  inclusion  of  the  station  delay  from  the  Pn  inversion  provides  an 
additional  constraint  on  the  Moho  depth  and  average  crustal  P  velocity  in  the  joint 
inversion  with  receiver  functions  and  surface  wave  dispersions. 

3.  Technical  Approach 

We  have  implemented  a  joint  inversion  method  that  includes  surface  wave  dispersion 
data,  receiver  functions,  and  Pn  station  delay  time.  A  feature  of  this  joint  inversion  is  that 
all  the  individual  components  (surface  wave  dispersion,  receiver  function,  and  Pn  delay 
time)  are  quantified  with  respect  to  the  same  station  of  interest.  Below  we  describe  in 
detail  the  joint  inversion  method  (see  also  Xu  et  al.,  2013a). 

3.1  Sensitivity  tests 

The  combination  of  receiver  functions  and  surface  wave  dispersion  bridges  some 
resolution  gaps  associated  with  each  individual  data  set  (e.g.,  Julia  et  al.  2000;  Ozalaybey 
et  ah,  1997).  Fig.  1  shows  an  example  of  joint  inversion  of  receiver  functions  and  surface 
wave  dispersion  at  an  INDEPTH  station.  However,  important  ambiguity  remains.  In  the  P 
receiver  function,  the  Moho  P-to-S  converted  phase  (Ps)  is  generally  the  strongest.  The 
Moho  depth  determined  from  the  Ps  delay  (relative  to  P)  trades  off  strongly  with  the 
crustal  Vp/Vs  ratio.  Multiple  converted  phases  would  help  resolving  the  ambiguity  (Zhu 
and  Kanamori,  2000),  but  they  are  usually  more  difficult  to  identify.  Sensitivity  tests 
(Fig.  2)  suggest  that  adding  well-constrained  P  information  will  help  the  joint  inversion 
of  receiver  function  and  surface  wave  dispersion  to  find  the  right  S  model  when  the 
Vp/Vs  ratio  is  not  known,  which  is  often  the  case.  One  type  of  P  information  is  the  Pn 
station  delay  time,  which  provides  constraint  on  Moho  depth  and  average  crustal  velocity. 
This  is  the  type  of  P  information  we  focus  on  in  this  project. 

3.2  Implementation  of  the  joint  inversion  method 

We  have  successfully  implemented  a  joint  inversion  method  that  includes  surface  wave 
dispersion  data,  receiver  functions,  and  Pn  delay  times  that  are  associated  with  the  same 
station  of  interest.  The  implementation  includes  the  following  key  elements.  (1)  We  use 
fast  global  search  method,  the  neighborhood  algorithm  (NA)  (Sambridge,  1999).  Using  a 
global  parameter  search  approach,  compared  with  a  more  traditional  linearized  iterative 
inversion  approach,  has  two  key  advantages:  (a)  We  have  more  control  on 
parameterization,  for  example,  to  include  Pn  data,  we  need  to  define  a  Moho;  and  (b)  we 
are  able  to  explore  parameter  space  more  fully  and  to  quantity  the  resolution  and 
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uncertainty  of  the  solution  more  easily.  (2)  Our  parameterization  includes  the  Moho,  the 
depth  of  which  is  allowed  to  vary.  It  can  include  variable  crustal  layers  as  well.  Regions 
between  layers  are  parameterized  as  smooth  functions  (linear  or  splines)  to  reduce  the 
number  of  search  parameters.  (3)  The  implementation  uses  parallel  programming  with 
MPI  calls,  which  greatly  speeds  up  the  model  search  with  multi-processor  or  parallel 
machines.  Fig.  3  shows  a  synthetic  example  using  the  global  search  method  for  joint 
inversion  of  surface  wave  and  receiver  functions  and  Pn  delay  time.  The  method  finds  a 
solution  with  perfect  fit  to  the  input  model  and  data. 

3.3  Data  preparation 

We  use  the  HiClimb  data  to  test  our  joint  inversion  implementation.  This  is  a  large 
volume  of  data  (some  180  stations),  and  we  need  a  effective  procedure  to  obtain  each 
data  set  for  the  joint  inversion.  First,  we  obtained  high  resolution  dispersion  curves  along 
the  array  using  a  linear  array  inversion  procedure.  Because,  it’s  a  near-linear  dense  array, 
we  can  reduce  the  3D  structure  inversion  problem  to  a  2D  problem  along  the  array,  thus 
allowing  us  to  refine  the  grid  size  (to  0.1  degree).  We  implemented  a  version  of  ambient 
noise  surface  wave  tomography  that  is  suitable  for  a  dense  linear  array.  Second,  we 
constructed  receiver  functions  for  all  the  stations  using  the  semi-automatic  procedure 
described  in  section  4.2.  The  interactive  procedure  allowed  us  to  control  data  quality  but 
also  to  increase  efficiency,  which  is  essential  to  processing  a  large  volume  of  data.  For 
each  receiver  function,  we  also  extract  the  delay  time  of  the  Moho-converted  phase  Ps 
relative  to  the  reference  phase  P  from  the  H-K  stacking  procedure  (Zhu  and  Kanamori, 
2000).  Although  the  crustal  thickness  from  the  receiver  function  trades  off  strongly  with 
the  Vp/Vs  ratio,  the  Ps  delay  time  is  robust  as  long  as  the  converted  phase  is  strong 
enough.  The  delay  time  provides  an  efficient  constraint  on  the  model  search  by  ruling  out 
models  with  significant  differences  in  the  delay  time.  Third,  our  joint  inversion  method 
includes  Pn  station  delay  times  obtained  from  Pn  tomography  of  the  whole  of  western 
China  (see  below). 

3.4  Practical  considerations  for  an  implementation  for  real  data 

Below  we  describe  practical  considerations  in  the  implementation  for  real  data  (with 
HiClimb  data  as  an  example).  Particular  effort  was  made  in  formulating  effective 
procedures  and  in  model  parameterization.  Although  a  search-based  approach  provides 
great  flexibility  in  fitting  different  types  of  data,  it  is  generally  very  expensive 
computationally.  The  selection  of  parameterization  requires  particular  care  because  of  the 
need  to  limit  the  number  of  parameters. 

Derivation  of  starting  model.  Our  joint  inversion  is  done  in  two  steps.  In  the  first  step, 
we  invert  only  the  dispersion  data  to  obtain  a  starting  S  velocity  profile  at  each  station, 
which  is  then  used  as  the  starting  model  for  the  model  searches.  The  main  purpose  of  the 
starting  model  is  to  reduce  the  model  space  to  help  the  search  algorithm  converge  faster. 
We  use  cubic  splines  for  model  parameterization.  The  advantage  of  using  splines  is  that 
they  yield  a  naturally  smooth  model  and  yet  reduce  the  number  of  parameters  for 
searching. 
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For  the  dispersion  inversion  in  the  first  step,  we  adopt  12  uniformly  distributed  spline 
nodes  from  the  surface  down  to  150  km.  We  seek  the  solution  that  satisfies  the  following 
minimizer: 

min  (||Dobs-Dpred||+/l||m||+  <p||Lm||)  (1), 

where  the  first  term  is  data  misfit,  Dobs  is  the  data  (dispersion  curve)  we  are  trying  to  fit, 
and  Dpred  is  the  predictions  from  the  random  model  generated  by  the  NA  algorithm.  The 
second  term  is  the  model  constraint.  Since  the  dispersion  periods  of  the  HiClimb  data 
from  the  ambient  noise  extends  only  up  to  45  s,  which  has  limited  constraint  on  the 
deeper  part  of  the  structure,  we  introduce  a  model  constraint  for  the  depth  from  95  km 
down  to  150  km  using  an  S  model  of  China  and  the  surrounding  region,  which  includes 
group  velocity  period  up  to  120  s  (Xu  et  al.,  2013b).  The  third  term  is  the  smoothing 
constraint.  L  is  the  Laplacian  smoothing  operator  (Lees  and  Crosson,  1989),  which 
requires  the  second  derivative  of  the  model  to  be  zero.  X  and  <p  are  two  real  positive 
numbers  that  have  been  chosen  so  that  different  terms  are  well  balanced.  We  seek  the 
solution  that  minimizes  the  LI  norm  of  above  formula  instead  of  the  traditional  L2  norm. 
The  choice  of  the  LI  norm  reduces  the  risk  that  the  model  is  biased  by  extreme  values 
due  to  error  in  the  data.  We  convert  the  spline  model  into  a  layered  velocity  model  to 
calculate  predicted  Rayleigh  group  and  phase  velocities.  We  divide  the  whole  150  km 
into  30  layers  with  5  km  thickness  for  each  layer.  Velocity  within  each  layer  is  constant. 
We  fix  the  Vp/Vs  ratio  for  all  the  layers  to  1.75.  Density  is  calculated  by  the  Nafe-Drake 
(1963)  relationship.  The  forward  calculation  subroutine  is  adopted  from  Robert 
Herrmann’s  program  (1991). 

Joint  inversion.  The  parameterization  for  the  second  step  in  the  joint  inversion  is 
different  from  the  first  step.  We  use  two  sets  of  cubic  splines  to  represent  the  velocity 
model:  one  for  the  crust  (first  set)  and  one  for  the  Moho  and  the  mantle  (second  set).  The 
last  node  of  the  first  spline  has  the  same  value  as  the  first  node  of  the  second  spline  so 
that  the  model  is  continuous  throughout  the  depth.  The  node  spacing  for  each  set  of  spline 
is  not  uniform;  instead,  it  is  stretched  over  depth,  with  closer  spacing  between  nodes  at 
the  top  of  the  spline  and  larger  spacing  at  the  bottom  of  the  spline.  We  use  the  following 
formula  to  map  the  normalized  spline  space  (0  to  1)  to  the  depth: 

h  =  Dx2+h0  (2), 

where  h  is  the  depth  from  the  top  of  the  first  node  in  the  spline;  D  is  the  total  thickness;  ho 
is  the  depth  of  the  first  node;  x  is  the  coordinate  in  the  normalized  spline  space  between  0 
and  1.  This  set  up  is  effective  in  modeling  the  structure  at  shallow  depth  and  near  the 
Moho  discontinuity  region  by  giving  more  freedom  to  the  model  where  it  is  better 
constrained  by  short  period  dispersion  data  and  receiver  functions,  respectively.  It  also 
yields  a  smoother  model  at  the  lower  crust  and  deeper  mantle,  which  are  less  constraint 
by  the  data.  We  use  12  nodes  for  the  first  spline  in  the  crust  and  10  nodes  for  the  second 
spline  in  the  Moho  and  the  mantle.  We  also  allow  the  depth  of  the  connecting  node 
between  the  first  and  the  second  splines  to  vary  so  that  the  Moho  depth  can  be  adjusted. 
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The  use  of  the  Moho  and  mantle  spline  provides  more  flexibility  compared  with  a  Moho 
with  an  imposed  discontinuity  or  linear  gradient.  For  a  Moho  depth  at  70  km,  the  node 
spacing  between  the  first  and  second  nodes  in  the  mantle  is  -0.99  km  and  the  spacing 
between  the  second  and  third  nodes  is  -2.96  km,  enough  to  model  a  sharp  velocity  jump, 
thus  mimicking  the  Moho  discontinuity.  On  the  other  hand,  this  parameterization  also 
works  well  in  the  situation  where  the  Moho  is  not  sharp  but  presents  as  a  transition  zone. 
We  convert  the  spline  model  into  layered  model  as  in  the  first  step.  We  use  2  km  layer 
thicknesses  throughout  the  depth  range  of  the  first  spline  in  the  crust.  In  the  mantle,  we 
use  fifteen  2  km  layers  for  the  top  30  km  as  a  transition  zone  followed  by  5  km  layers  to 
the  bottom  at  150  km.  We  seek  the  solution  that  satisfies  the  LI  norm  of  the  following 
minimizer: 


min(w!  IlDSl’sp-D^p  ||  +  w2 1| D",bss  -  D|?frsed ||  +  w3  ||D^fay  -  D||a®ady||  +  A||m  -  m0||  +  cp||Lm||) 

(3), 

where  the  first  term  is  the  dispersion  misfit,  the  second  term  is  the  RFs  misfit  and  the 
third  and  last  terms  are  the  same  as  formula  (1).  The  wi,  W2  and  W3  are  the  relative 
weighting  for  dispersion,  RFs  and  Pn  delay  time  misfits,  respectively. 

We  impose  model-screening  algorithms  to  speed  up  model  convergence.  Every  time  NA 
generates  a  random  model,  we  compare  the  model  with  the  starting  model  from  the 
dispersion  inversion  for  each  layer  except  those  15  layers  in  the  transition  zone.  If  the 
velocity  difference  in  any  given  layer  is  larger  than  0.5  km/s,  this  model  is  discarded  by 
setting  a  large  misfit  value  without  going  through  the  forward  calculation.  We  also 
require  the  velocity  to  increase  monotonically  from  the  last  two  nodes  in  the  crust  to  first 
two  nodes  in  the  mantle.  This  model  screening  process  can  reduce  computational  time  by 
confining  the  model  search  in  the  model  space  where  the  likelihood  of  the  model  is  the 
largest. 

Another  screening  parameter  is  the  measured  Ps  delay  time.  Because  the  converted  Ps 
phase  may  be  mis-identified  in  the  H-K  stacking  process,  we  need  to  sort  out  the 
“outliers”  of  the  Ps  delay  times  before  applying  the  screening.  We  fit  the  measured  Ps 
delay  times  of  all  the  stations  with  a  smooth  spline  curve.  If  the  measurement  for  a  station 
is  different  from  the  prediction  of  the  spline  fit  by  more  than  1.5  s,  then  we  use  the  spline 
prediction  as  the  Ps  delay  time  for  that  station,  otherwise  we  use  the  original  Ps  delay 
measurement.  For  each  model  generated  in  the  search,  we  calculate  the  velocity  gradient 
from  each  layer  to  the  next  layer  and  pick  the  one  with  largest  gradient  as  the  “Moho”.  In 
most  cases,  it  is  the  last  layer  in  the  first  spline  to  the  first  layer  in  the  mantle  spline.  We 
calculate  the  predicted  Ps  delay  time  using  the  “Moho”  defined  above  and  compare  it 
with  the  measured  Ps  delay.  If  the  difference  is  larger  than  1  s,  then  the  model  is 
discarded  (assigned  a  large  misfit  value).  Our  experience  suggests  that  by  introducing  the 
Ps  delay  time,  the  model  converges  faster  than  without  the  Ps  delay  screening. 

The  density  of  each  layer  is  calculated  using  the  same  Nafe-Drake  relationship  as 
described  above.  Because  of  the  lack  of  coherent  crust  multiple  signals,  only  the  first  15  s 
of  receiver  functions  are  included  in  inversion.  For  each  station,  we  select  4  receiver 
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functions  with  different  ray  parameters  for  inversion.  The  move  out  of  Ps  signal  at 
different  ray  parameters  is  helpful  in  determining  the  Moho  depth.  Subroutines  of 
dispersion  and  receiver  function  forward  calculation  are  adopted  from  Robert  Herrmann’s 
computer  programs.  We  iterate  800  times  to  obtain  the  final  model  and  100  new  models 
are  generated  each  iteration  in  5  Voronoi  cells. 

3.5  Computation  considerations 

The  project  involves  a  tremendous  of  amount  of  data  and  computation,  so  special 
considerations  on  computational  efficiency  and  computational  resources  are  needed. 

First,  we  spent  considerable  effort  in  speeding  up  the  calculations  of  empirical  Green 
functions  of  surface  waves  from  noise  correlations.  Because  the  task  for  each  pair  is  quite 
independent,  it  can  be  parallelized  efficiently.  The  main  obstacle  is  lots  of  I/O  operations, 
which  require  expensive  reading  and  writing  hard  drives.  We  have  come  up  a  solution 
that  would  speed  up  the  computation  by  distributing  tasks  to  multi-processors  within  the 
same  computer  (node)  and  execute  them  simultaneously.  The  method  uses  the  LINUX 
POSIX  message  queue  (mq)  capabilities.  Second,  we  developed  a  semi-automatic 
procedure  to  generate  receiver  functions.  Because  of  large  number  of  stations  and  large 
amount  of  data,  such  a  procedure  is  essential  for  productivity  (see  below).  Third,  the  joint 
inversion  is  based  on  MPI-calls,  which  can  be  run  on  parallel  machines.  We  used  a  Linux 
cluster  in  our  institution  and  the  National  Center  for  Supercomputing  Applications  on  the 
University  of  Rhode  Island  campus  for  most  of  the  computations. 


4.  Results  and  Discussion 

We  now  summarize  briefly  the  results  of  the  individual  components  of  surface  wave 
dispersion,  receiver  functions,  and  Pn  tomography  and  describe  in  some  detail  the  joint 
inversion  results. 

4.1  Dispersion  measurements,  dispersion  maps,  and  3D  S  model 

We  have  systematically  measured  dispersion  for  the  construction  of  dispersion  maps 
from  ambient  noise  correlations  as  well  as  traditional  earthquake  data.  The  stations 
include  the  Chinese  national  backbone  stations  and  open  permanent  and  temporary 
stations  from  the  early  1990s  with  a  total  of  over  500  stations  (Fig.  4).  Earthquake  data 
include  events  greater  than  magnitude  5.0  from  2000  to  2007.  The  total  dataset  has 
dispersion  measurements  for  the  period  from  8  s  to  120  s  with  a  maximum  of  -35000 
measurements  for  group  velocity  and  -22500  measurements  for  phase  velocity  at  the 
period  of  20  s. 
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Dispersion  maps  at  different  periods  are  constructed  by  minimizing  the  following 
equation: 

min  (||Am— r|p+A^||m|p+  <^||Lm|p) 

The  first  term  is  data  misfit,  in  which  m  is  the  parameter  vector,  A  is  the  coefficient 
matrix,  r  is  the  vector  of  travel  time  residuals.  The  second  term  is  the  model 
regularization  with  a  damping  parameter  of  the  non-zero  real  number  A.  The  third  term  is 
the  smoothing  constraint  with  the  non-zero  real  number  cp.  The  smoothing  operator  L  is 
the  Laplacian  operator  (Lees  and  Crosson,  1989). 


Examples  of  new  dispersion  maps  are  shown  in  Fig.  5.  Earthquake  data  account  for  about 
50%  of  the  total  group  velocity  measurements  for  periods  from  8  to  70  s.  We  compare 
group  velocity  dispersion  maps  from  earthquake  data  and  the  results  using  the  whole  data 
set.  All  large-scale  major  features  are  similar,  e.g.  low  velocity  in  sedimentary  basins  at 
short  periods,  a  low  velocity  zone  throughout  the  Tibetan  Plateau,  and  an  east- west 
velocity  contrast  at  mid  periods,  although  there  are  differences  in  detailed  structures. 


We  invert  the  dispersion  maps  for  a  3D  S  velocity  model.  Figs.  6-7  show  some  S  velocity 
maps  and  profiles,  respectively.  A  detailed  discussion  of  the  main  features  of  the  model 
can  be  found  in  Xu  et  al.  (2013b). 


4.2  Receiver  Functions 

We  have  systematically  processed  receiver  functions  of  stations  in  Western  China  (Fig. 
8),  using  a  semi-automatic  procedure  we  have  developed  to  generate  receiver  functions. 
The  procedure  includes  interactive  review  and  selection  of  good  quality  raw  data  using  a 
graphic  interface,  picking  and  alignment  of  P  waveforms  using  waveform  cross¬ 
correlation  and  a  graphic  interface,  automatic  calculations  of  receiver  functions  and 
stackings  according  to  ray  parameters  and  azimuths,  and  interactive  selections  of  good 
quality  receiver  functions  using  a  graphic  interface. 

4.3  Pn  Tomography 

We  have  worked  systematically  on  Pn  waves.  Our  data  come  from  three  sources:  1) 
Chinese  bulletin  data,  2)  international  bulletin  data  (ISC),  and  3)  hand-picked  data  at 
temporary  and  permanent  stations.  Distribution  of  permanent  stations  in  western  China 
(even  including  regional  and  provincial  networks  in  China)  is  quite  sparse.  However, 
there  have  been  a  number  of  PASSCAL  experiments  in  Western  China.  These  stations 
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are  usually  in  remote  areas  where  permanent  stations  are  lacking,  providing 
complimentary  ray  coverage. 

For  hand-picking,  we  designed  an  inter-active  scheme  to  pick  Pn  arrival  times  semi- 
automatically.  It  involves  inspecting  original  seismograms  (vertical  component,  filtered 
at  0.2  to  3  Hz),  discarding  noisy  traces,  zooming  in  Pn  arrivals,  picking  each  trace  by 
hand,  plotting  travel-time  curves  and  reduced  travel-time  curves,  and  double-checking 
picks  if  there  are  obvious  anomalies.  We  have  examined  about  200  shallow  events  (depth 
less  than  40  km)  with  magnitude  greater  than  5.5  over  15  years  (1995-2009)  in  the  region 
from  25  to  56  deg  north  and  68  to  108  deg  east.  We  require  epicentral  distances  to  be 
from  2  to  12  degrees  and  at  least  10  Pn  picks  for  each  event.  We  picked  over  1800  Pn 
arrivals  from  81  events. 

We  combine  the  three  data  sources  and  impose  the  following  data  selection  criteria  for 
the  combined  data  set.  (1)  Earthquakes  whose  depth  is  over  45  km  are  discarded.  (2)  The 
apparent  velocity  of  rays  is  limited  between  7.0  and  9.0  km/s.  (3)  To  ensure  robust 
constraint  on  the  station  delays,  we  keep  the  stations  with  at  least  5  Pn  picks  and  keep  the 
earthquakes  with  at  least  5  Pn  records.  (4)  We  use  an  earthquake-grouping  technique 
(Liang  et  al.,  2004)  to  locate  earthquake  clusters,  then  choose  the  earthquake  with  the 
largest  number  of  Pn  picks  in  a  cluster  and  exclude  all  the  other  earthquakes  in  the  same 
cluster.  After  the  data  selection,  we  obtain  54,849  rays  from  398  stations  and  6,220 
earthquakes. 

For  Pn  tomography,  we  follow  the  basic  approach  of  parameterizing  the  Pn  residual  into 
mantle  perturbation,  the  earthquake  term,  and  the  station  term  (e.g.,  Hearn,  1996;  Liang 
et  al.,  2004;  Liang  and  Song,  2006).  For  the  initial  model  of  inversion,  we  choose  average 
crustal  velocity  of  6.3  km/s,  average  Pn  velocity  of  8.06  km/s,  and  average  Moho  depth 
of  51  km  (Liang  and  Song,  2006).  We  divide  the  study  region  into  1°  x  1°  2-D  grids. 

The  inversion  results  of  Pn  velocities  and  station  delays  are  shown  in  Fig.  9.  The  Pn 
velocity  map  (Fig.  9 A)  shows  a  number  of  significant  features  consistent  with  the 
geology  of  the  region  and  previous  studies,  including  high  Pn  velocity  anomalies  in  four 
major  basins  (Sichuan,  Tarim,  Qaidam,  and  Junggar  basin),  high  uppermost  mantle 
velocities  in  the  southern  Tibetan  Plateau  (TP),  low  Pn  velocities  in  northern  TP,  and  a 
striking  low  velocity  belt  extending  from  west  to  east  beneath  northern  TP  to  the  eastern 
and  southeastern  margins  of  the  plateau,  observed  previously  by  Liang  and  Song  (2006). 
The  station  delays  vary  greatly  in  the  study  region  (Fig.  9B),  suggesting  large  changes  in 
crust  thickness,  as  we  expect. 
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The  station  delays  provide  a  critical  data  set  for  our  joint  inversion  with  surface  waves 
and  receiver  functions.  We  can  examine  the  constraint  of  the  Pn  station  delay  time  by 
examining  the  crustal  thickness  estimated  using  the  station  delay  alone. 


The  station  delay  is  related  to  the  crustal  velocity  ( vc ),  crustal  thickness  ( H ),  and  the 
reference  model  (vc,  H ) 

V/2 


as 


follows:  t„—t- -  =  H 

'i  n 


/  Vc  /  Vm 


^2  “  /-2 
Vc  /  V„ 


where  77  = 


1/  _  1/ 

7  2  7  2 

7  V  /  v 


,  is  the  vertical  slowness  in  the  crust.  If  we  assume  an  average 

V/  vc  /  vmJ 

crustal  velocity,  we  can  estimate  the  crustal  thickness.  Fig.10  shows  a  comparison  of  the 
crustal  thickness  along  HiClimb  stations  estimated  from  three  data  sets:  Pn  station  delay, 
receiver  function,  and  joint  inversion  of  surface  wave  dispersion  and  receiver  function. 
We  see  that  the  crustal  thickness  derived  from  different  data  sets  is  mutually  consistent  in 
general.  There  are,  however,  some  discrepancies,  which  likely  arise  from  the 
uncertainties  in  the  crustal  velocities  and  which  we  hope  to  resolve  with  joint  inversions 
of  all  the  data  sets. 


4.4  Joint  Inversions 


Combining  individual  components  described  above,  we  now  focus  on  joint  inversions  of 
surface  wave  dispersion,  receiver  function,  and  Pn  delay  time. 

1)  Tests  on  joint  inversions. 

We  first  conduct  various  tests  of  the  joint  inversions.  Key  considerations  include:  (1)  the 
inclusion  of  the  Vp/Vs  ratio  in  the  search  parameter  in  the  coding,  (2)  the  relative 
weighting  between  different  types  of  data  (dispersion,  receiver  function,  and  Pn  delay 
time),  (3)  the  difference  of  model  results  with  or  without  the  Pn  data,  (4)  the  robustness 
of  the  joint  inversion. 

Fig.  11  shows  an  example  of  the  joint  inversion  for  a  HiClimb  station.  The  joint  inversion 
with  Pn  shows  a  very  similar  S  velocity  profile  to  that  of  the  joint  inversion  with  surface 
wave  dispersion  and  receiver  function  only.  The  predictions  fit  to  the  observed  dispersion 
data  and  to  receiver  functions  very  well.  The  fits  are  not  distinguishable  from  the  results 
of  the  joint  inversion  using  the  dispersions  and  receiver  functions  only.  However,  we  note 
that  there  is  a  noticeable  difference  in  Moho  depth  due  to  a  difference  in  the  best-fit 
Vp/Vs  ratio  (1.83  compared  to  a  assumed  value  of  1.75  in  the  dispersion  and  receiver 
function  inversion).  The  example  underlines  the  importance  of  including  P  constraints  in 
the  joint  inversion. 


A  certain  question  to  ask  is  how  robust  the  derived  Vp/Vs  ratio  from  the  joint  inversion 
is.  We  thus  conduct  another  test  with  different  weightings  on  the  Pn  data  in  the  total 
misfit  function  (Fig.  12).  Our  tests  show  that  the  derived  Vp/Vs  ratio  is  quite  robust  for  a 
range  of  Pn  data  weightings  when  the  quality  of  receiver  functions  is  good.  However,  the 
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sharpness  of  the  Moho  varies.  Essentially,  the  Vp/Vs  ratio  is  constrained  by  the  timing  of 
dispersion,  the  Ps  delay  times  in  the  receiver  functions,  and  the  Pn  delay  time.  However, 
the  sharpness  of  the  Moho  affects  the  amplitudes  of  the  Ps  in  the  receiver  functions.  A 
greater  Pn  weighting  down-weighs  the  receiver  function  in  the  misfit,  thus  the  receiver 
functions  are  not  fit  as  well. 


Yet  another  test  is  the  effect  of  the  Vp/Vs  ratio  and  the  inclusion  and  exclusion  of  Pn  data 
on  the  model  inversion,  shown  in  Fig.  13  for  a  HiClimb  station  HI 350.  In  previous 
inversion  with  the  availability  of  the  Pn  delay  time,  we  set  the  initial  Vp/V s  value  to  be 
1.75  and  the  inversion  result  from  surface  wave  and  receiver  functions  is  shown  in  blue. 
With  the  addition  of  the  Pn  delay  time,  we  found  the  optimal  Vp/Vs  value  of  1.88.  As  a 
result,  the  inversions  (red  and  green)  with  the  optimal  Vp/Vs  of  1.88  yield  significantly 
different  Moho  depth  and  uppermost  mantle  structures  than  the  initial  inversion. 

However,  the  inclusion  of  Pn  delay  time  also  down-weights  the  dispersion  and  receiver 
function  data  in  the  total  misfit  function.  As  a  result,  the  joint  inversion  with  the  Pn  data 
included  yields  a  more  gradual  Moho  (green)  and  smaller  Ps  amplitude  in  the  receiver 
function  (Fig.  13).  However,  if  we  use  the  new  Vp/Vs  ratio  of  1.88  and  include  only  the 
dispersion  and  receiver  function  data,  we  obtain  a  similar  S  velocity  profile  as  in  the  joint 
inversion  as  well  as  a  sharper  Moho,  which  fits  the  receiver  function  Ps  amplitude  better. 

2)  Strategies  for  the  joint  inversions 

From  the  above  tests,  we  outline  our  strategies  for  the  joint  inversion  of  real  data.  For  a 
“good”  station  that  has  all  the  three  types  of  data  (surface  wave,  receiver  function,  and  Pn 
delay  time)  and  the  Moho  conversion  is  strong  in  the  receiver  functions,  we  first  use  all 
the  data  to  derive  an  optimal  Vp/Vs  ratio,  and  then  fix  the  Vp/Vs  ratio  at  the  optimal 
value  and  conduct  a  separate  inversion  using  dispersions  and  receiver  functions  only.  For 
other  stations,  we  use  our  best  estimates  of  the  Vp/Vs  ratios  based  on  the  “good”  stations 
and  perform  the  joint  inversions  of  dispersions  and  receiver  functions  only. 

3)  Joint  inversion  results:  HiClimb  stations 

Fig.  14  shows  Vp/Vs  ratios  derived  from  HiClimb  stations.  The  results  that  use  all  the 
stations  that  have  the  three  types  of  data  (about  60  stations)  (open  circles)  are  consistent 
with  those  from  a  subset  of  stations  that  have  the  “best”  receiver  functions  (with  strong 
Moho  conversions)  (solid  dots).  The  derived  Vp/Vs  ratios  show  large  variations  from  the 
collision  front  across  the  Tibetan  Plateau  (TP).  The  Vp/Vs  ratios  in  the  TP  is  significantly 
greater  than  in  the  Indian  Shield.  However,  the  ratios  in  the  Lhasa  Block  are  significantly 
smaller  than  in  the  Himalaya  and  Qiangtang  Blocks. 

Fig.  15  shows  results  of  the  joint  inversion  of  the  three  types  of  data  for  HiClimb 
stations.  The  joint  inversion  (Fig.  15  A)  is  generally  consistent  with  the  joint  inversion  of 
dispersions  and  receiver  functions  only  (Fig.  15B)  or  the  inversion  of  the  dispersions 
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only  (Fig.  15C).  But  there  are  some  noticeable  differences  in  mid-crustal  low  velocity 
zone  structure  and  in  Moho  depth. 


4)  Joint  inversion  results:  western  China 

Using  the  strategy  outlined  above,  we  have  systematically  performed  joint  inversions  for 
each  station  that  has  all  the  three  types  of  data  (surface  wave,  receiver  function,  and  Pn 
delay  time).  From  the  joint  inversion,  we  obtain  an  optimal  Vp/Vs  ratio  for  each  station. 
These  Vp/Vs  ratios  are  then  interpolated  to  derive  a  Vp/Vs  map  (Fig.  16).  Note  because 
of  non-uniform  distribution  of  stations  for  which  we  can  perform  the  joint  inversions,  the 
derived  Vp/Vs  ratio  map  is  highly  uncertain.  However,  we  do  observe  consistent  large- 
scale  patterns.  In  particular,  the  Vp/Vs  ratios  are  anomalously  high  for  western  TP 
compared  to  surrounding  regions. 

An  important  limitation  of  the  joint  inversion  method  is  that  it’s  limited  by  station 
distributions.  However,  because  of  excellent  coverage  of  surface  wave  dispersions  for  the 
3D  S  velocity  model  and  our  derived  Vp/Vs  map,  we  can  derive  a  3D  P  velocity  model 
(Fig.  17)  based  on  the  3D  S  velocity  model  and  the  Vp/Vs  map.  The  derived  3D  P  model 
shows  distinct  geological  features  of  tectonic  blocks,  which  mimics  features  from  the  3D 
S  model  with  some  modifications.  In  the  future,  it  would  be  useful  to  compare  and 
integrate  such  a  derived  P  model  with  P  travel-time  derived  models  for  the  region. 


5.  Conclusions 

We  successfully  developed  and  implemented  a  search-based  algorithm  (NA)  for  joint 
inversion  of  surface  wave  dispersion  data,  receiver  functions,  and  Pn  delay  times  using 
parallel  programming.  We  obtained  unprecedented  data  sets  for  surface  wave  dispersion, 
receiver  functions,  and  Pn  tomography,  needed  for  the  joint  inversion.  We  tested  the  joint 
inversion  method  and  proposed  a  practical  strategy  for  joint  inversions  of  the  real  data. 
We  obtained  detailed  joint  inversion  results  for  the  dense  HiClimb  array,  which  are 
generally  consistent  with  previous  results  but  show  considerable  difference  in  details 
(crustal  structure  and  Moho  depths).  We  have  systematically  performed  joint  inversions 
of  all  available  stations  in  western  China  and  obtained  a  3D  S  velocity  model,  an  average 
crustal  Vp/Vs  map,  and  a  derived  3D  P  velocity  model  of  the  region. 
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Fig.  1.  An  example  of  commonly-used  joint 
inversion  of  receiver  functions  and  surface  wave 
dispersion.  The  left  panel  shows  the  S  velocity 
model  (black  lines)  from  the  inversion  and  the  H- 
k  stacking  result  (red  lines).  The  right  panel 
shows  receiver  function  and  surface  wave 
dispersion  fits,  black  represents  data  and  red 
represents  prediction  from  the  obtained  S 
velocity  model.  The  Vp/Vs  ratios  are  fixed  using 
the  values  obtained  from  the  H-k  stacking 
results. 


A.  Assume  input  average  Vp  (6.25  km/s) 


B.  Assume  wrong  Vp  (6.5  km/s) 


Vs  (km/s)  t  (s) 


Vs  (km/s)  t  (s) 


Fig.  2.  Synthetic  tests  on  combining  Pn  station  delays  with  receiver  functions  and  surface  wave  dispersion 
for  a  joint  inversion.  The  input  model  is  listed  in  Table  1.  (A)  Correct  average  crustal  Vp  (6.25  km/s)  is 
assumed,  which  would  yield  the  correct  Moho  depth  (30  km)  from  Pn  station  delays.  (B)  Average  crust  Vp 
is  assumed  to  be  a  wrong  value  of  6.5  km/s,  yielding  Moho  depth  of  33.3  km  from  Pn  station  delays.  The  S 
structure  recovered  in  B  is  considerably  poorer  than  in  A  and  has  a  Moho  (slightly  less  than  30  km) 
inconsistent  with  the  Pn  Moho  (33.3  km). 
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(a) 


S  velocity  (km/s) 


Fig.  3.  Synthetic  test  of  joint  inversion  including  receiver  functions,  surface  wave  dispersion  and  Pn  delay 
time  data  using  the  Neighborhood  Algorithm.  The  input  velocity  model  is  shown  in  blue  in  (a).  The  input 
model  consists  of  three  layers.  The  Moho  discontinuity  is  at  30  km  depth.  The  Vp/Vs  ratios  in  the  crustal 
layers  are  1.75  and  1.85,  respectively.  Two  sets  of  splines  are  used  in  the  inversion.  Each  set  consists  of  5 
evenly  spaced  nodes.  The  output  velocity  model  is  plotted  in  red  in  (a).  Note  the  velocity  jump  in  the  crust  is 
smoother  than  in  the  input  model  due  to  the  natural  smoothness  of  the  spline  function.  An  average  Vp/Vs 
ratio  of  1.80  in  the  crust  is  obtained  during  searches  when  Pn  delay  time  data  are  included.  The  input 
(solid  line)  and  output  (dashed  line)  of  receiver  function  waveforms  are  shown  in  (b).  Only  the  first  10 
seconds  of  the  receiver  functions  are  used  in  the  inversion.  Rayleigh  wave  group  and  phase  dispersion  data 
(blue)  from  5  s  to  45  s  are  used  in  the  inversion  (c).  The  corresponding  dispersion  predictions  from  the 
output  model  are  shown  in  red. 
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Fig.  4.  Distributions  of  seismic  stations  (triangles)  and  earthquakes  (stars)  used  for  joint  ambient  noise  and 
earthquake  surface  wave  tomography. 


Approved  for  public  release;  distribution  is  unlimited. 

13 


50 


40" 


3C 


20' 


2.7 


I 

2.8  2.9  30 

Group  velocity  (km/a) 

(a) 


27  2  0  2,9  3  0  3-1  3.2 

Group  velocity  (km/s) 

(b) 


BO"  90"  IDO'  110"  120'  130" 


80  90'  1O0"  110"  120'  13d' 


3.4 

3  5  3-6  3  7 

3-B 

34 

3S 

3  6  3.7 

3.8 

Group  velocity  (km/s) 

Group  velocity  (km/s) 

(e) 

(f) 

Fig.  5.  Group  velocity  maps  at  10,  30,  and  50  s  using  EGFs  and  earthquake  data  (left)  and  earthquake  data 
only  (right).  Results  from  two  different  data  sets  are  generally  consistent. 
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Fig.  6.  S  velocity  maps  at  different  depths  from  our  inverted  model.  White  great  circle  paths  in  (a)  indicate 
the  cross  sections  shown  in  Fig.  7. 
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Fig.  7.  S  velocity  cross-sections  of  our  model  in  western 
China.  The  panels  (a)  through  (g)  correspond  to  the  cross 
sections  along  great  circle  paths  labeled  in  Fig.  6a. 
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Fig.  8.  Stations  for  which  receiver 
functions  have  been  systematically 
processed  and  obtained  for  the  joint 
inversion  of  this  project. 
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Fig.  9.  Results  of  preliminary  Pn  tomography  in  western  China.  (A)  Pn  velocity  map.  (B)  Station  delays. 
The  number  of  stations  is  around  400,  including  bulletin  stations  and  PASSCAL  stations. 
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Fig.  10.  Crustal  thickness  along  the  profile  of  HI-CLIMB 
stations.  Circles  are  Moho  depths  calculated  from  the  Pn 
station  delays  for  each  station.  The  red  line  is  the 
approximated  Moho  location  from  the  receiver  function 
study  by  Nabelek  et  al.,  (2009)  using  HI-CLIMB 
stations.  The  blue  line  is  the  approximated  Moho 
location  from  joint  inversion  of  receiver  function  and 
surface  wave  dispersion  by  Xu  et  al.  (2013).  The 
topography  of  Tibetan  Plateau  is  shown  on  the  top  of 
the  figure. 
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Fig.  11.  Test  result  of  joint  inversion  with  Pn  data  on  HiClimb  station  HI  620.  a)  S  velocity  profile  result 
with  Pn  data  (red  line  with  the  obtained  Vp/Vs  and  the  fit  to  Pn  delay  time  labeled),  compared  with  the 
result  without  Pn  (blue  line);  b)  The  fit  for  surface  wave  dispersion  curves  (observation  in  blue  and 
prediction  in  red),  including  group  velocity  (lower  values)  and  phase  velocity  (higher  values);  c)  The  fit  for 
4  receiver  functions  with  different  ray  parameters  (receiver  function  data  are  in  solid  lines  and  predicted 
waveforms  are  in  dashed  lines). 
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Fig.  12.  Tests  of  Pn  weighting  on  the  inversion  results  (using  HiClimb  station  HI  620).  The  blue  curve  is 
the  result  from  inversion  of  dispersion  and  receiver  functions.  The  red  curves  are  joint  inversions  with  the 
addition  of  Pn  data  with  different  Pn  weightings  as  labeled. 
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Fig.  13.  Tests  of  model  sensitivities  to  Vp/Vs  ratio  and  with/without  Pn  data  constraint.  The  example  is  for 
HiClimb  station  HI 3 50.  The  initial  Vp/Vs  value  is  set  to  1. 75.  The  optimal  value  of  Vp/Vs  from  surface 
wave  dispersion,  receiver  functions,  and  Pn  station  delays  is  found  to  be  1.88. 
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Fig.  14.  Vp/Vs  ratios  derived  from  joint  inversion  of  surface 
wave  dispersion,  receiver  functions,  and  Pn  station  delays. 
The  stations  are  from  the  HiClimb  array.  The  curves  are 
spline  fits  to  the  data.  Solid  dots  show  the  results  from  the 
stations  with  the  “best”  Moho  conversions  (a  total  of  32 
stations).  Open  circles  show  the  results  for  all  the  stations 
that  have  all  three  types  of  data  (a  total  of  60  stations). 


A) 


B) 


3.5  4.0  4.5 


Vs  (km/s) 


E 

£ 

Q_ 

<d 

Q 


-200  -100  0  100  200  300  400  500 


-0.3 


-0.2  -0.1  -0.0  0.1  0.2  0.3 

Amplitude 


Fig.  15.  Results  of  joint  inversion  for  the  Hi-CLIMB  array  (A,  B).  Crosses  and  circles  indicate  the  Moho 
depth  under  the  stations.  Crosses  indicate  a  strong  Moho  discontinuity  while  circles  mark  a  gradual  Moho 
velocity  transition.  (A)  shows  the  joint  inversion  with  the  three  types  of  data  and  (B)  shows  the  joint 
inversion  with  dispersion  and  receiver  functions  only.  Shown  in  comparison  are  results  from  surface  wave 
dispersions  only  (C)  and  from  common  conversion  point  stacking  of  the  receiver  functions  only  (D).  The 
thick  white  line  is  the  approximated  Moho  location  from  the  receiver  function  study  by  Nabelek  et  al. 
(2009)  using  the  same  Hi-CLIMB  stations. 
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Fig.  16.  Map  of  Vp/Vs  ratios  in  western  China 
obtained  in  this  study.  The  symbols  are  from  joint 
inversions  and  colors  are  from  interpolation. 


Fig.  17.  Maps  of  the  P  velocity  model  of  western  China  from  this  project. 
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